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Abstract 



This paper presents a stability test for a class of interconnected nonlinear systems motivated 
by biochemical reaction networks. The main result determines global asymptotic stability of 
the network from the diagonal stability of a dissipativity matrix which incorporates information 
about the passivity properties of the subsystems, the interconnection structure of the network, 
and the signs of the interconnection terms. This stability test encompasses the secant criterion 
for cyclic networks presented in [1], and extends it to a general interconnection structure rep- 
resented by a graph. The new stability test is illustrated on a mitogen activated protein kinase 
(MAPK) cascade model, and on a branched interconnection structure motivated by metabolic 
networks. The next problem addressed is the robustness of stability in the presence of diffusion 
terms. The authors use a compartmental model to represent the localization of the reactions 
and present conditions under which stability is preserved despite the diffusion terms between 
the compartments. 

1 Introduction 

This paper continues the development of passivity-based stability criteria for interconnected sys- 
tems motivated by classes of biochemical reaction networks. In [1,2] the authors studied a cyclic 
interconnection structure in which the first subsystem of a cascade is driven by a negative feedback 
from the last subsystem downstream. This cyclic feedback structure is ubiquitous in gene regulation 
networks [3-14], cellular signaling pathways [15, 16], and has also been noted in metabolic path- 
ways [17,18]. In [1,2] the authors first presented a passivity interpretation of the "secant criterion" 
developed earlier in [8, 14] for the stability of linear cyclic systems, and next used this passivity 
insight to extend the secant criterion to nonlinear systems. The notion of passivity evolved from 
an abstraction of energy conservation and dissipation in electrical and mechanical systems [19,20], 
into a fundamental tool routinely used for nonlinear system design and analysis [21,22]. 

The first contribution of this paper is to expand the analysis tool of [1] to a general interconnec- 
tion structure, thus obtaining a broadly applicable stability criterion that encompasses the secant 
criterion for cyclic systems as a special case. As in [1], our approach is to exploit the passivity 



properties and the corresponding storage functions [19] for smaller components that comprise the 
network, and to construct a composite Lyapunov function for the interconnection using these stor- 
age functions. The idea of using composite Lyapunov functions has been explored extensively in 
the literature of large-scale systems as surveyed in [23,24], and led to several network small-gain 
criteria [25, 26] that restrict the strength of the interconnection terms. A distinguishing feature 
of our passivity-based criterion, however, is that we take advantage of the sign properties of the 
interconnection terms to obtain less conservative stability conditions than the small-gain approach. 

To determine the stability of the resulting network of passive subsystems we follow the formalism 
of [27,28], and construct a dissipativity matrix (denoted by E below) that incorporates information 
about the passivity properties of the subsystems, the interconnection structure of the network, and 
the signs of the interconnection terms. As a stability test for the interconnected system, we check 
the diagonal stability [29] of this dissipativity matrix, that is, the existence of a diagonal solution 
D > to the Lyapunov equation E T D + DE < which, if feasible, proves that the network is 
indeed stable. In particular, the diagonal entries of D serve as the weights of the storage functions 
in our composite Lyapunov function. Although similar results can be proven by combining the pure 
input /output approach in [27,28] with appropriate detectability and controllability conditions (see 
the discussion in Section [H] below) , the direct Lyapunov approach employed in this paper allows us 
to formulate verifiable state-space conditions that guarantee the desired passivity properties for the 
subsystems. These conditions are particularly suitable for systems of biological interest because 
they are applicable to models with nonnegative state variables, and do not rely on the knowledge 
of the location of the equilibrium. 

The second contribution of this paper is to accommodate state products which are disallowed 
in the nonlinear model studied in [1]. This is achieved with a new storage function construction 
for each subsystem which, in the absence of state products, coincides with the construction in [1]. 
Thanks to this extension, our stability criterion is now applicable to a broader class of models, even 
in the case of cyclic systems. This class encompasses a mitogen activated protein kinase (MAPK) 
cascade model with inhibitory feedback proposed in [15,16], which is studied in Example [Q as an 
illustration of our main result. The final result in the paper employs a compartmental model to 
describe the spatial localization of the reactions, and proves that, if the passivity-based stability 
criterion holds for each compartment and if the storage functions satisfy an additional convexity 
property, then stability is preserved in the presence of diffusion terms between the compartments. 

The paper is organized as follows: Section [2] gives an overview of the main results in [1]. Section 
[3] presents a general interconnection structure represented by a graph, and gives the main stability 
result of the paper. Section U] illustrates this result on biologically motivated examples. Section O 
studies robustness of stability in the presence of diffusion terms in a compartmental model. Section 
[6] develops an extension of the purely input /output-based proof in [2] of the secant criterion to the 
general graphs studied in this paper. In doing so, it adapts a lemma from [27, 28] and compares 
the input/output approach with the state-space results derived earlier in the paper. Section [7] gives 
the conclusions. 
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2 Overview of the Secant Criterion for Cyclic Systems 



To evaluate stability properties of negative feedback cyclic systems, references [8, 14] analyzed the 
Jacobian linearization at the equilibrium, which is of the form 
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1, • • • ,n, and showed that A is Hurwitz if the following sufficient condition 
b\ ■ ■ ■ b n 



< sec(7r/n) r 



(2) 



Unlike a small-gain condition which would restrict the right-hand side of ([2]) to be 1, the "secant 
criterion" ([2]) also exploits the phase of the loop and allows the right-hand side to be as high as 8 
(when n = 3). The secant criterion is also necessary for stability when the a^s are identical. 

Local stability of the equilibrium proven in [8, 14], however, does not rule out the possibility 
of periodic orbits. Indeed, the Poincare-Bendixson Theorem of Mallet-Paret and Smith for cyclic 
systems [30,31] allows such periodic orbits to coexist with stable equilibria, as we illustrate on the 
system: 



X3 



-XX + cp(x 3 ) 

-X2 + Xx 
-x% + X 2 



(3) 



where 

tp(x 3 ) = e- 10 ^-!) + .lsat(25(x 3 - 1)), (4) 

and sat(-) := sgn(-) min{l, |-|} is a saturation^ function. The function is decreasing, and its slope 
has magnitude 63 = 7.5 at the equilibrium xx = x 2 = x 3 = 1. With ax = a 2 = a 3 = bx = b 2 = 1 
and n = 3, the secant criterion ([2]) is satisfied and, thus, the equilibrium is asymptotically stable. 
However, simulations in Figure [T] show the existence of a periodic orbit in addition to this stable 
equilibrium. 

To study global stability properties of cyclic systems with negative feedback, in [1,2] the authors 
first developed a passivity interpretation of the secant criterion ([2]), and next used this passivity 



One can easily modify this example to make <p(-) smooth while retaining the same stability properties. 
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Figure 1: Trajectory of (J3j) starting from initial condition x = [1.2 1.2 1.2] T , projected onto the 
x\-X2 plane. 



insight to extend the secant criterion to the nonlinear model: 

x\ = -fi(xi) + h n {x n ) 
±2 = ~f2(x 2 ) + hi(x{) 

Xn — fn{Xn) ~\~ h n — \ \X n — l) 



(5) 



in which Xi G JR>o, fi(-), i = 1, • • ■ , n and /ij(-), i = 1, • • • , n — 1 are increasing functions, and fo n (-) 
is a decreasing function which represents the inhibition of the formation of x\ by the end product 
x n . When an equilibrium x* exists, [1] proves its global asymptotic stability under the following 
condition: 



dhi(xj) 
dxi 



dxi 

7i • • - 7„ < sec(7r/n) n , 



>o> 



1, ■ • ■ , n, 



(6) 
(7) 



which encompasses the linear secant criterion ([2]) with ji = bi/ai. 

A crucial ingredient in the global asymptotic stability proof of [1] is the observation that the 
secant condition is necessary and sufficient for the diagonal stability of the matrix 



E, 



cyclic 



-1/71 
1 -1/72 
1 







-1/73 





■l/7n 
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that is, for the existence of a diagonal matrix D > such that 



E. 



•T 

cyclic 



D + DE cydlc < 0. 



(9) 
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The authors of [1] connect this diagonal stability property to the global asymptotic stability of 
^ by first noting that the assumption ([BJ) guarantees an output strict passivity property, where 
l/7i quantifies the excess of passivity in each Xj-subsystem. They then incorporate this passivity 
information in the diagonal terms of the dissipativity matrix ([8]) , and represents the interconnection 
structure with the off-diagonal terms. Finally, they use the diagonal stability condition ([9]), which 
is equivalent to the secant criterion ([7]), to check whether the excess of passivity in each subsystem 
overcomes the loss of passivity in the interconnection. In particular, the diagonal entries of D 
constitute the weights of the storage functions in a composite Lyapunov function for ([5]) . 

3 From the Cyclic Structure to General Graphs 

We now extend the diagonal stability procedure outlined above for cyclic systems to a general 
interconnection structure, described by a directed graph without self-loops. The nodes represent 
subsystems with possibly vector outputs, and a separate link is used for each output channel. For 
the nodes i = 1, • ■ ■ ,N and links I = 1, ■ ■ ■ , M, we denote by Cf C {1, • • • , M} the subset of links 
for which node i is the sink, and by the subset of links for which node i is the source. We 
write i = source(Z) if I £ , and i = sink(Z) if I £ Cf. Using this graph we introduce the dynamic 
system: 



where xi £ M>o, and /«(•), <7i(-)> * = 1) "** ,N, hi(-), I = 1, ••• ,M are locally Lipschitz functions 
further restricted by the following assumptions: 

Al: fi(0) = and, for all a > 0, g^a) > 0, hi(a) > 0. 

Assumption Al guarantees that the nonnegative orthant iR> is forward invariant for (|10p . The 
strict positivity of gi(xi) is also essential for our analysis since we exploit the sign properties of 
^K^sourcc^)) which are multiplied by g%(xi) in ([10]) . 

A2: There exists an equilibrium x* £ -ZR>q for (jlOp . 

A3: For each node i, the function fi(xi) / g%(xi) satisfies the sector property: 



A4-' For each node i, and for each link I £ £■ , the function hi(xi) satisfies one of the following 
sector properties for all Xi G M>o — {x*}: 



To distinguish between positive and negative feedback signals we assign to each link I a positive 
sign if (fT2|) holds, and a negative sign if (fT3|) holds, and rewrite (fT2|) - (fT3|) as 




(10) 




(11) 



fa -x*)[hi{ Xi ) -hi(x*)} > 
(xi -x*)[hi(xi) -hi(x*)} < 0. 



(12) 
(13) 



sign(link l)(x{ - x*)[hi(xi) - h(x*)] > 



(14) 
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Vxi £ M> - {x*}. 

A5: For each link I E C~ there exists a constant 7/ > such that, Mxi E iR>o — {%*}, 

■ tv 1 n h l( x i) ~ h l{x*) . ( . 

S1 gn(hnk Q — _ m) < 7l . (15) 

Theorem 1 Consider the system (jlOp . and suppose assumptions A1-A5 hold. If the M x M 
dissipativity matrix 

r -i/ 7 j if k=i 

i?^ = < sign(link A;) if source(Z) = sink(/c) (16) 
^ otherwise 

is diagonally stable; that is, if there exists a diagonal matrix D > such that 

E T D + DE < 0, (17) 

then the equilibrium x* is asymptotically stable. If, further, for each node i one of the following 
two conditions holds, then x* is globally asymptotically stable in M> : 
a) C~ is nonempty and there exists at least one link I £ CJ such that 



px 

lim / 



fl'i(o-) 

b) is empty; that is, the outdegree of node i is zero 



Xi hAa) - hi(x*) , x 

v l - = 00, (18) 



'^t4 = oo, (19) 



X,'— >oo 



and there exists a class- /Coo function™] such that 



{x l -x\){^- f -^^>\x i -x\\^{\x i -x\\) y Xl >0. (20) 
\9i\Xi) 9i\ x i)/ 

Proof: We first prove the theorem for the case when is nonempty for all i = 1, • • ■ ,iV; that 
is, when there are no nodes with outdegree equal to zero. In this case we construct a composite 
Lyapunov function of the form 

M 

V(x - X *) = J2 diVi(x sourcc(l) - < ourcc(0 ) (21) 
1=1 

in which the components are 



-source! 



e(i) hi(a) - hi(x somce{l) ) 
m^sourceCO ~ ^sourcef/)) = sign(link I) / — da (22) 



2 fC is the class of functions M>o — > -R>o which are zero at zero, strictly increasing and continuous. KLoc is the 
subset of class-/C functions that are unbounded. 
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and the coefficients di > are to be determined. The function (|2ip is positive definite because each 
component VJ is a positive definite function of (x sc , U rce(0 — x l OUTCe ^) due to the sign property (fT3j) of 
the integrand in ([22]) . and because (x SOUICC ^ — x* outcc ^) = 0, I = 1, • • ■ , M, guarantees x — x* = 
by virtue of the fact that each node is the source for at least one link. 

We now claim that the function V\ in ([22|) satisfies the dissipativity property 

M 

Vi<yiY, E ikVk (23) 
k=l 

where 

yi := Sign(link [^source©) ~ M<ource(i))] ( 24 ) 

I = 1, • • ■ , M, and the coefficients i^fc are as in (I16p . Before we prove this claim, we first note 
that the diagonal stability property (fT7|) and the estimate (f23|) together imply that the Lyapunov 
function (|2ip . with coefficients d[ obtained from the diagonal elements of D, yields a negative definite 
derivative from which asymptotic stability of x* follows. If, further, for each node i there exists at 
least one link Z G C~ such that (fl~8l) holds, then the Lyapunov function (j2~T1) grows unbounded as 
\x\ — > oo, thus proving global asymptotic stability. 

If there exist nodes with outdegree equal to zero, then the arguments above prove that the 
subsystem comprising of the nodes with outdegree one or more is asymptotically stable. The out- 
puts hi from this subsystem serve as inputs to the nodes with outdegree equal to zero. Because 
the dynamics of these nodes are of the form (|10p and are asymptotically stable by A3, asymptotic 
stability for the equilibrium x* follows from standard results on cascade interconnections of asymp- 
totically stable systems (see e.g. [32, p. 275]). Likewise, when condition (b) holds, (|20p and (|19p 
imply an input-to- state stability (ISS) property [33] for the driven subsystem of the cascade, and 
global asymptotic stability follows because the cascade interconnection of an ISS system driven by 
a globally asymptotically stable system is globally asymptotically stable [33]. 

We conclude the proof by showing that the claim (f23|) is indeed true. To this end we compute 
from (|22p and (|10p the derivative 

Vi = sign(lmk l^hixi) - h{x*)] f-lipl + U A (25) 
m := ^ h k {x sourcc{k) ). (26) 



where i = source (Z), and 



Adding and subtracting 



k£C 



within the bracketed term in (|25|) . we obtain 



_n(p) +u ._ * 

gi(xi) gi{x*) 
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Next, noting that sign(link l)[hi(x{) — hi(x*)] and ( l^'l ~ TTnAl P ossess the same signs due to 



IT]) and (I14p . and using (I15p . we obtain the inequality 



- sign(hnk 0[/»,(xi) - /i,«)] f ^4 " 44rl ^ "-fafa) " HO?- (29) 

U^i) 9i{ x *)J 7* 

Substituting ([59]) in (|28p . and using the variables y\ defined in ([23]) . we get 

^ = --i/? + yi(«t-0- (30) 

7/ 

Finally, noting from ([26]) and ([27]) that 

Ui - u* = sign(link fe)^, (31) 

we rewrite ()30p as 

Vj < + y« V sign(link k)y k , (32) 

fce£+ 

which is equivalent to (|23p by the definition of the coefficients E^i in (j!6p . □ 

The assumptions A3-A5 rely on the knowledge of the equilibrium x* which may not be available 
in practice. When the functions fi(-), 9i(-), and hi(-) are C 1 , the following incremental conditions 
guarantee A3-A5, and do not depend on x*: 

A3': For each i = !,-■■ ,N, and Vxj > 0, 



d ( fi(xi) 



dxi \gi(xi) 
A4 ': For each / = !,-■■ ,M, and Va;* > 0, 



> 0. (33) 



signflink > 0. (34) 

OXi 

A5': For each link I € there exists a constant 72 > such that 

<li Vxi>0. (35) 





9hi(xi) 






dxi 






( fi(xi)\ 





Although the growth assumption (|35[) may appear restrictive, most biologically relevant nonlin- 
earities satisfy this condition globally. If there exist closed intervals Xi C M>o such that X\ x ■ ■ ■ x^jy 
is forward invariant for (jlOj) . a less conservative 7; may be obtained by evaluating (|35p on A^, rather 
than for all Xj > 0. This relaxation is particularly useful in biological applications where Xi repre- 
sents the amount of a substance which may be lower- and upper-bounded. 

The dissipativity matrix E in (j!6p combines information about the interconnection structure of 
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the network with the passivity properties of its components. Because the off-diagonal components 
of this matrix are negative for links that represent inhibitory reaction rates, diagonal stability is less 
restrictive than a networked small-gain condition [25,26] which ignores the signs of the off-diagonal 
terms. In the case of a cyclic graph where each link I = 1, • • • , n connects source i = I to sink 
i = l + l (modn), and where only link n has a negative sign, (|16p assumes the form ([8]). Theorem Q] 
thus recovers the result of [1] as a special case, and further relaxes it by accommodating the gi(xi) 
functions in (|1U|) which are not allowed in [1]. 



4 Examples 

Example 1 To illustrate Theorem [1] we first study a simplified model of mitogen activated protein 
kinase (MAPK) cascades with inhibitory feedback, proposed in [15,16]: 

di(l-xi) n 

ci + xi ei + (1 — xi) 1 + kx 3 

b 2 x 2 d 2 (l-x 2 ) , , 

x 2 = + —j- rX! (37) 

c 2 + x 2 e 2 + (1 - x 2 ) 

b3X3 d 3 (l-x 3 ) 

x 3 = — + , n ^x 2 . (38) 

C3 + X3 e 3 + (1 - x 3 ) 

The variables Xi E [0, 1] denote the active forms of the proteins, and the terms 1 — xi indicate 
the inactive forms (after nondimensionalization and assuming that the total concentration of each 
of the proteins is 1). The second term in each equation indicates the rate at which the inactive 
form of the protein is being converted to active form, while the first term models the inactivation 
of the respective protein. For the proteins Xi, i = 2,3, the activation rate is proportional to the 
concentration of the active form of the protein upstream, which facilitates the conversion. 
The activation of the first protein x±, however, is inhibited by X3 as represented by the decreasing 
function fJ./(l + kx 3 ). 

The model d36D-([3HD is of the form ([TOD with 



fi(Xi) — 1 gi{Xi) — . , i — 1, 2, 3, 

Q T" Xi Ci ~r yl Xi ) 

hi{xi) =x h i = 1,2, h 3 {x 3 ) = 1 _^ ■ (39) 

Because the underlying graph is cyclic with each link / = 1,2,3 connecting source i = I to sink 
i = l + l(mod3), and because /i3(-) is strictly decreasing, the dissipativity matrix E in (|16p is of 
the form and, as proved in [1], its diagonal stability is equivalent to the secant criterion ([7]). 
However, unlike the model © of [1] which disallows state products, Theorem Q] above accommodates 
the functions gi(xi), and is applicable to ([36]) - ([38]) . 

To reduce conservatism in the estimates for the 7j's in Theorem[T]we further restrict the intervals 
[0, 1] in which x^s evolve by noting that ^3(3:3) takes values within the interval Because 
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^3(^3) is the input to the x\ -subsystem, and because the function 6i : [0, 1] — > [0, 00) defined by 

fiixi) 



i(xi) ■-- 



(40) 



is strictly increasing, it follows from the bounds on the input signal that the interval X\ = 
[xi t min, xi jma x] '■= + k)),9^ 1 (fi)] is an invariant and attractive set for the xi-subsystem. 

Since x\ and X2 serve as inputs to the x 2 - and ^-subsystems respectively, the same conclusion 
holds for the intervals X 2 = [x2,min,x 2 ,max] and X 3 = [x 3tmin ,x 3tmax ], where 



'i,max 



i—l,max ) 



e 3 = 0.01, 



(41) 



Xi.min ■ — $j {Xi— l,mm) 

i = 2, 3. With the following coefficients from [34]: 

61 = ei = ci = b 2 = 0.1, c 2 = e 2 = c 3 
63 = 0.5, ^1 = ^2 = ^3 = 1, [i = 0.3, 



we obtained 7i's numerically by maximizing the left-hand side of (|35[) on Xi for various values of the 
parameter k. This numerical experiment showed that the secant condition 717273 < 8 is satisfied 
in the range k < 4.35 (for k = 4.36 we get 717273 = 11.03). Reference [34] gives a small-gain 
estimate k < 3.9 for stability, and shows that a Hopf bifurcation occurs at around k = 5.1. The 
estimate k < 4.35 obtained from Theorem Q] thus reduces the gap between the unstable range and 
the small-gain estimate. 

Example 2 A common form of feedback inhibition in metabolic networks occurs when several 
end metabolites in different branches of a pathway inhibit a reaction located before the branch 
point [18,35]. As an example of this situation we consider the network in Figure [2] where the end 
metabolites with concentrations #4 and xq inhibit the formation of x\ from an initial substrate xq. 
Assuming that xq is kept constant, and that its conversion to x\ is regulated by two isofunctional 



x\ 



X'2 




X3 



xq 



Figure 2: Feedback inhibition in a branched network. The dashed links 4 and 7 indicate negative 
(inhibitory) feedback signals. The dissipativity matrix obtained from (|16f) for this network is (14311 . 
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enzymes each of which is selectively sensitive to £4 or xq, we represent this network as in fjlOp : 



x\ 

■h 

X3 
X4 

X5 
X6 



-fi(xi) + h4,(x 4 ) + h 7 (x 6 ) 
-hfa) + h\{x{) 
-hfa) + h 2 (x 2 ) 
-U(xi) + ^3(2:3) 
-f5(x 5 ) + 715(2:2) 
-feixe) + h(x 5 ), 



(42) 



where the functions ^4(^4) and hj(xe) are decreasing due to the inhibitory effect of X4 and xq, 
while hi(-), I = 1, 2, 3, 5, 6 and /*(•)> i = 1, ■ ■ ■ ,6 are increasing. 

Rather than study specific forms for these functions, we assume that Al and A2 hold, and that 
7z's exist as in (|35p . An application of Theorem [1] then proves global asymptotic stability of the 
equilibrium if the dissipativity matrix 
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_ J_ 

76 
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_ J_ 



(43) 



is diagonally stable. Note that the 4x4 principal submatrices obtained by deleting row-column 
pairs {5,6,7} and {2,3,4} each exhibit a cyclic structure for which, as shown in [1], diagonal 
stability is equivalent to the secant criteria 



71727374 < sec(7r/4) 4 = 4 and 71757677 < 4, 



(44) 



respectively. Because principal submatrices of a diagonally stable matrix are also diagonally stable, 
we conclude that (|44|) is a necessary condition for the diagonal stability of (|43|). In fact, we prove 
the following necessary and sufficient condition: 



Lemma 2 The matrix E in |^3[ ) is diagonally stable iff 

71727374 + 71757677 < sec(7r/4) 4 = 4 . 



(45) 



Proof: We prove the sufficiency of this condition as a consequence of a more general fact. Consider 
the following diagonal matrix: 



D = diag 1 , 



7374 74 



72 7273 



7677 77 2 \ 
2 ' 75 ' 7576/ 



(46) 
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and the matrix 



M := E T D + DE. 



We will prove that condition (|45[) implies that M < 0. Diagonal stability of E follows from this 
claim in view of the following argument: Given any 7j's satisfying the constraint (|45p . we can 
find 7i > 7i that still satisfy the constraint, and under this transformation E gets transformed to 
£7 = E + A, where A is some positive diagonal matrix. Now let -D be defined for E as in (|46p 
with 7j's replaced by 7j's. Since £ T Zi> + DE < E T D + DE = M, and since Af < 0, it follows that 
E T D + DE <0, which means that E is diagonally stable. 

To prove that (|45p implies Af < 0, we let E e := E — el for each e > 0, and show that 
M e = EjD + D-E e is negative definite for small enough e > 0. By continuity, this last property 
implies that M < 0. In order to check negative definiteness of M e , we consider the principal minors 
//i(e), i = 1, . . . , 7 of M e , and ask that they all have sign (—1)* for small e > 0. Each m is a 
polynomial of degree < 7 on e. The determinant of M e can be expanded as follows: 

(£) = 87477(75 + 2 76 + 77X72 + 2 7 3 + 74) Ag2 + 0(e3) (4?) 
7172737576 

where A = 71727374 + 7i7s7677 — 4. Similarly, we have: 

M £ ) = " 274772( V 2 ? +74) A g + 0( g 2 ), 
71727375 

(£) = 2747677(72 + 273 + 74) Ag + ^ 
71727375 

/ \ -274(72 + 273 + 74) a e , n / 2N 

iu, 4 {£) = 5 Aie + (J(e ), 

717273 



where Ai = 71727374 — 4, 



7l 



and 



M3(e) = T i± ^Ai + 0{e), 
27172 

^ 2 ( e) = Z232i(Ai-4) + O(e), 
47172 



2 

= 2e. 

7i 



Since Ai < A, we conclude that the matrix M £ is negative definite for all small enough e > if 
and only if A < 0. In particular, condition (|45p implies that M < 0, as claimed. 

Finally, we prove the necessity of (|45|) for the diagonal stability of E in (|43|) . To this end, we 
define 2£ = diag (71, • • • ,77) E which has all diagonal components equal to —1, and characteristic 
polynomial equal to: 

( a + l) 3 [( S + l) 4 + fc], 

where k := 71727374 + 71 757677- For k > 0, the roots of (s + l) 4 = —k have real part ±y/k/4 — 1; 
hence A; < 4 is necessary for these real parts to be negative. Because (|45p is necessary for the 
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Hurwitz property of E, it is also necessary for its diagonal stability. Since diagonal stability of E is 
equivalent to diagonal stability of E, we conclude that (|45|) is necessary for the diagonal stability 
of E. 

5 Stability of a Compartmental Model with Diffusion 

A compartmental model is appropriate for describing the spatial localization of processes when each 
of a finite set of spatial domains ("compartments") is well-mixed, and can be described by ordinary 
differential equations. Instead of the lumped model (|1U|) . we now consider n compartments, and 
represent their interconnection structure with a new graph in which the links k = 1, • • • , m indicate 
the presence of diffusion between the compartments j = 1, • • • , n they interconnect. Although the 
graph is undirected, for notational convenience we assign an orientation to each link and define the 
n x m incidence matrix S as 

{+1 if node j is the sink of link k 
— 1 if node j is the source of link k (48) 
otherwise. 

The particular choice of the orientation does not change the derivations below. 

We first prove a general stability result (Theorem [3] below) for a class of compartmental models 
interconnected as described by the incidence matrix S. We then apply this result in Corollary 2] to 
the situation where the individual compartments possess dynamics of the form studied in Section 
El We let 

Xj := (a^i, • • • ,Xj jN ) T 

be the state vector of concentrations Xj^ in compartment j, and let Xj = Fj(Xj) represent the 
dynamics of the jth compartment in the absence of diffusion terms. Next, for each link k = 
1, • • • , m, we denote by 

fJ>k,i( x 

sink(fc),i •^source(fc),i ) (49) 

the diffusion term for the species i, flowing from source(A;) to sink(A;), and assume the functions 
Mfc,i(')> k = 1, ■ •■ ,m, i = 1, • •• ,N, satisfy 

vVk,i( a ) ^ °> Vo- G M. (50) 

Then, the coupled dynamics of the compartments become: 

Xj = Fj{Xj) + (Sj, I N )M(S T ® I N )X) j = 1, • • • , n (51) 

where Sj,. is the jth row of the incidence matrix S, In is the N x N identity matrix, "£§>" represents 
the Kronecker product, 

X:=[X?...XZ\ T (52) 
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and n : R mN -> R mN is denned as 

H(z) := [/ii,iOi) • • • Hi,n(zn) 



Mm,l( z (m-l)7V+l) ' " " ^m,N{z m N)Y. 



(53) 



We now prove stability of the coupled system (j5Tj) under the assumption that a common Lya- 
punov function exists for the decoupled models Xj = Fj(Xj), j = 1, • • • , n, and that this common 
Lyapunov function consists of a sum of convex functions of individual state variables: 

Theorem 3 Consider the system (|5ip where the function //(•) is as in (|53|) and (|50p . If there exists 
a Lyapunov function V : R N — > M of the form 

F(x) = y 1 (xi) + --- + tV(ziv) (54) 

where each Vi(sj) is a convex, differentiable and positive definite function, satisfying 

VV(x)F j (x) < -a(\x\) j = l,---,n (55) 

for some class-/C function a(-), then the origin X = of (|5ip is asymptotically stable. If, further, 
V(-) is radially unbounded, then X = is globally asymptotically stable. 



Proof: We employ the composite Lyapunov function 

n 

V{X) = ^2v(Xj), 

i=i 

and obtain from (|5TT) and (1551): 



(56) 



V(X) < - ct(\Xj\) + [W(Xi) • • • W(X n )](S ® / 7V ) / u((5 T ® /at)X). 
We next rewrite the second term in the right-hand side of (|57p as 



and note from (|4"8]) that ([58]) equals 



W T (Xi) 



w T pr n ) 



T 



fi((S T ®I N )X), 



source (fc) ; 



fc=l 



(57) 



(58) 



(59) 



where (ik,it i = 1, • • ■ , AT, denotes the diffusion function (|49f) . and the argument is dropped for 
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brevity. Next, using ([33]) . we rewrite ([59]) as 



m N 

i)] A*fc,i- (60) 

k=l i=l 

Because Vj(-) is a convex function, its derivative VV^(-) is a nondecreasing function and, hence, 
Wi( •^sink(fc),j a; sourcc(fc),i) possesses the same sign as (£ s ink(fc),i ~~ ^source(fc),i)- We next recall 

from the sector property (f50l) that the function /i/^j in (fl9l) possesses the opposite sign of its 
argument (x s ink(fc),i ~~ x sourcc(fc),i)- This means that each term in the sum ([60|) is nonpositive and, 
hence, ([57]) becomes 

n 

V(aO<-X>(l*il), (61) 

from which the conclusions of the theorem follow. □ 

Theorem [3] is applicable when each compartment is as described in Section [3l hi(-) satisfies (|34p . 
and <7i(-)'s, i = 1, • • • , N, are nonincreasing functions. This is because the Lyapunov construction 
(|2ip in Section [3] consists of a sum of terms as in (|54p , each of which is convex when the derivative 
of (|22p is nondecreasing: 

Corollary 4 Consider the system (|51|) where the function /i(-) is as in (|53|) and (|5U|) . and Fj(x), 
j = 1, • ■ ■ ,n, are identical and represent the right-hand side of (jlOp . If all assumption of Theorem 
CD hold and if, in addition, hi(-) satisfies ([34"]) . and <7i(-)'s, i = 1, ■ • • , N, are nonincreasing functions, 
then the equilibrium X = [x* T , ■ ■ ■ -,x* T ] T is globally asymptotically stable. 

6 Comparison of the State-Space and Input/Output Approaches 

The earlier paper [2] gave a purely input/output (instead of state-space) version of the secant 
criterion, phrased in the language of passivity of L 2 operators. We now explain how to extend this 
I/O approach to the general graphs studied in this paper. The result follows easily by imposing an 
appropriate diagonal stability condition, combined with a key lemma due to Moylan and Hill [27], 
and Vidyasagar [28]. Below we give a streamlined version of this lemma, and compare it with the 
state space approach employed earlier in this paper. 

We denote by L\ the extended space of signals (thought of as time functions) w : [0, oo) — > R 
which have the property that each restriction wt = H[o,T] is m 

L 2 (0,T), for every T > 0. Given 
an element w G L 2 and any fixed T > 0, we write \\w\\ T for the L 2 the norm of this restriction wt, 
and given two functions v,w € L 2 and any fixed T > 0, the inner product of vt and wt is denoted 
by (v,w)t- The same notation is used for vector functions. 

We view the M subsystems to be interconnected as operators Ej : L 2 — > L 2 : m i— > yi, and 
impose the following strict passivity property: there exist constants ji > ("secant gains" in [2]) 
such that 

WUiWr — li{yii u i)T for each i = 1, . . . , M and each T > . (62) 
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We then consider the interconnection where 



Ui(t) = Vi (t) + A iy (t) , (63) 

or just u = v + Av, where the v^s are external inputs, y = col(yi, . . . ,um), v = col(i>i, . . . ,vm), 
and the Ai, i = 1, . . . ,M are the rows of an interconnection matrix A G M. MxM . In other words, 
the ith subsystem receives as inputs an external input plus an appropriate linear combination of 
outputs from the remaining systems (including possibly feedback from itself, if the corresponding 
diagonal entry of A is nonzero). We introduce: 

E := A - r 

where V = diag (i, . . . , 

Lemma. Suppose that there exists a diagonal positive definite matrix D G jj-MxAf such that 

DE + E'D < . 

Then, the system obtained from the systems £, using the interconnection matrix A is 1? stable as 
a system with input v and output y. More precisely, there is some constant p > such that, for 
any u,v,y G (L^) M such that (|62j) and (f63l) hold, necessarily \\y\\ T < p\\v\\ T for all T > (and 
therefore also ||y|| < p||^||, if v G (L 2 ) M ). 

Proof: We pick an a > such that DE + £"D < —2a/, and observe that, for any T > and any 
function z G L 2 (0,T), it holds that 

r-T f-T i f-T 

|2 



{Dz,Ez) = [ z{s)'DEz{s)ds = [ -z' (s)(DE+E' D)z(s) ds < -a [ z'(s)z(s)ds 
Jo Jo 2 j 



-a \\z\ 



Fix an arbitrary T > 0, and write D = diag (d\, . . . , djv/)- Since, for each i, (y^, Uj — ^-yi)T > 0, it 
follows that also [diy^Ui — ^-yi}r > 0, or, in vector form: 

(£>y,«-r y ) T >o. 

Substituting u = v + Ay, we obtain: (Dy, v + -E/y)r > 0, from which, using the Cauchy-Schwartz 
inequality: 

/^IMItIMIt ^ ( d V, v }t > -{Dy,Ey) T > a\\y\\ 2 T 

for some f3 > 0. So ||y|| T < p ||n|| T , with /) = ^, as desired. □ 

State-space stability results may be obtained as corollaries, by combining this I/O result with 
appropriate detectability and controllability conditions, as discussed in [2]. However, the direct 
Lyapunov approach employed earlier in this paper allowed us to formulate verifiable state-space 
conditions that guarantee the desired passivity properties for the subsystems. These conditions are 
particularly suitable for systems of biological interest because they are applicable to models with 
nonnegative state variables, and do not rely on the knowledge of the location of the equilibrium. 
The state-space approach further made it possible to prove robustness of our stability criterion in 
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the presence of diffusion terms. 



7 Conclusions 

We have presented a passivity-based stability criterion for a class of interconnected systems, which 
encompasses the secant criterion for cyclic systems [1] as a special case. Unlike the result in [1], 
we have further allowed the presence of state products in our model. Our main result (Theorem [1]) 
determines global asymptotic stability of the network from the diagonal stability of the dissipativity 
matrix (|16|) which incorporates information about the passivity properties of the subsystems, the 
interconnection structure of the network, and the signs of the interconnection terms. Although 
diagonal stability can be checked numerically with efficient linear matrix inequality (LMI) tools [36] , 
it is of interest to derive analytical conditions that make explicit the role of the reaction rate 
coefficients on stability properties. Indeed our earlier paper [1] showed that the diagonal stability 
of negative feedback cyclic systems is equivalent to the secant criterion of [8, 14]. In Example 2 we 
have derived a similar analytical condition for a branched cyclic interconnection structure. Further 
studies for deriving analytical conditions for practically important interconnection structures would 
be of great interest. Another research topic is to extend the stability result for compartmental 
models with diffusion in Section [5] to partial differential equation models. On this topic we have 
reported preliminary results applicable to cyclic systems in [37], and are currently studying more 
general interconnection structures. 
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